%This file aims to create IVs---averaging market share and market dominance
%indicators across all outside associated markets.

clear
tic

textFileName1=sprintf('sysyr.raw'); %the input file: the hsayr_id and sysyr_id and numhsa for involved multi-region system
textFileName2=sprintf('input4createIV.raw'); %the input file of basic info of all hospitals
sysyr=dlmread(textFileName1);  %the matrix: hsayr_id sysyr_id numhsa
H=dlmread(textFileName2);  %the matrix: sysyr_id ahayr_id hsayr_id prechs adj_vnd bdtot sysid nprofit profit teaching ncompet ratio perc_mcr perc_mcd
choice=1:13;
ns=numel(choice);


%%%%%%%%%%%%%%%%%%%%
%The first way is to find out the largest HSA according to how many
%member hospitals locate in.
exo_mktdom=zeros(size(sysyr,1),ns);
exo_mktsh=zeros(size(sysyr,1),ns);
for i=1:size(sysyr,1)
    tempsys=H(H(:,1)==sysyr(i,2),:);
    %%%%%%%%%%
    %The following is to compute the exogenous big hsa for a hospital
    %system
    syshsa=sort(tempsys(:,3));
    unique_syshsa=unique(syshsa);
    occur=histc(syshsa,unique_syshsa);
    bighsa=unique_syshsa(occur==max(occur));
    exo_bighsa=bighsa;
    %exo_bighsa(exo_bighsa(:)==sysyr(i,1))=[];
    if min(abs(exo_bighsa-sysyr(i,1)))==0
        exo_bighsa=[];
    end
  if isempty(exo_bighsa)==0 
      part_mktdom=zeros(numel(exo_bighsa),ns);
      part_mktsh=zeros(numel(exo_bighsa),ns);
      for m=1:numel(exo_bighsa)
          prechs=H(H(:,3)==exo_bighsa(m),4);
          prechs(prechs==1)=[];
          if isempty(prechs)==0
          occur=histc(prechs,choice);
          occur=reshape(occur,1,[]);           
          part_mktsh(m,:)=occur/sum(occur); %a vector 1 by 13
          temp_prechs=prechs;
          temp_prechs(prechs==2|prechs==13)=[];
          if isempty(temp_prechs)==0
             occur=histc(temp_prechs,choice)';
             mktdom=choice(occur==max(occur));
             Tchoice=choice';                     
             vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
             part_mktdom(m,:)=vmktff';            
          end 
          end
      end
      exo_mktdom(i,:)=mean(part_mktdom,1);
      exo_mktsh(i,:)=mean(part_mktsh,1);
  end 
end


iv_mktdom=[sysyr(:,1:2),exo_mktdom];
iv_mktdom(~any(iv_mktdom(:,3:end),2),:)=[];  %### UPDATED on 07/25/2019
iv_mktsh=[sysyr(:,1:2),exo_mktsh];
iv_mktsh(~any(iv_mktsh(:,3:end),2),:)=[];  %### UPDATED on 07/25/2019

iv_mktdom=sortrows(iv_mktdom,[1,2]);
iv_mktsh=sortrows(iv_mktsh,[1,2]);

hsayr=unique(sort(iv_mktdom(:,1)));

output_ivdom=zeros(numel(hsayr),ns);
output_ivmsh=zeros(numel(hsayr),ns);
for i=1:numel(hsayr)
    output_ivdom(i,:)=mean(iv_mktdom(iv_mktdom(:,1)==hsayr(i),3:end),1);
    output_ivmsh(i,:)=mean(iv_mktsh(iv_mktsh(:,1)==hsayr(i),3:end),1);   
end


textFileName11=sprintf('iv_dom_tmp.csv');
dlmwrite(textFileName11,[output_ivdom,hsayr],'delimiter', ',', 'precision', 8)
textFileName12=sprintf('iv_msh_tmp.csv');
dlmwrite(textFileName12,[output_ivmsh,hsayr],'delimiter', ',', 'precision', 8)







